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ABSTRACT 

In this paper we present the equations and some basic applications of relativistic SPH. 
The equations are generated under the assumption that there are no interaction terms 
between the fluid modelled and the background space-time, i.e. we are neglecting the 
perturbations to the metric produced by the modelled fluid. This corresponds to a 
stationary metric, which for this work, we further limit to be static as well. The 
equations use a new signal velocity term and artificial viscosity to smear out the 
effects of strong shocks and are tested here against Newtonian and relativistic shocks. 

Key words: methods: numerical - hydrodynamics ~ Shockwaves. 



1 INTRODUCTION 

The advent of x-ray telescopes with resolutions of Chandra 
and RXTE has opened the window into high energy as- 
trophysics. Although we now have increasingly complex and 
detailed data sets in these higher wavebands, the theory and 
our ability to understand these images is falling far behind. 
The reason for this is the complexities involved in solving 
the equations of relativistic fluid dynamics, an essential as- 
pect of high energy astrophysical phenomena. In general, 
the simulation of these systems requires the solution of the 
General Relativistic field equations as well as the solution of 
the equations of fluid dynamics. This is a complex and dif- 
ficult computationa l prob l em w hich has only been met with 
fimited success. (See iFonl i2003l) for a detailed review) 

However, there is a range of problems where the metric 
can be specified independently of the fiuid motion and it is 
to this class of problems that the present paper is directed. 
Of these problems the simplest involves motion in the flat 
space of special relativity and a typic al example is a rela- 
tivistic jet (see iMarti fc MiilleJ 1^9^). Even under these 
simplified conditions, these types of flow are notorious for 
the computational complexities that arise when one tries to 
model them. This is basically due to the strong shocks and 
the narrow, dense structures that form easily in the super- 
sonic flow. 

One approach to these problems is to set up a finite 
difference scheme and make use of the exact Riemann so- 
lutions for an ideal gas. This approach has been very suc- 
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cessful in the hands of iMarti fc MuUeJ (Il999ll and others 
( Schneider ct al. ( f 993) for example) who have obtained 
very accurate results for one and two dimensional shock 
wave problems involving shock tubes and jets. A disadvan- 
tage of these methods is that the Riemann solution must 
be known and for complex problems this may be difficult to 
obtain. For example, in a relativistic collision of two streams 
new particles may be created. The current relativistic Rie- 
mann solutions do not include such processes. Another ex- 
ample is the fiuid dynamical model of high energy nuclei 
collisions where the equation of state is a complicated func- 
tion of the densities and energies of the colliding fiuid. No 
relativistic Riemann solutions have been obtained for these 
equations of state, so that approximate Riemann solvers, of 
doubtful accuracy and stability must be used. 

An alternative approach is to avoid the discontinuities 
by smoothing them over some length scale, usually by an 
artificial viscosity term, and evolving the hydrodynamical 
equations numerically. If this is done on an Eulerian grid, the 
questions of resolving the fine structures mentioned above 
arise, as do the problems of the massive grid distortions 
which will be experienced in the vicinity of compact objects. 
These problems can be avoided by using Lagrangian parti- 
cle methods. The method which forms the basis of this pa- 
per, is the particle method SPH (for a review see Monaghaiil 
(|l992) and for applications to one dimensional shock prob- 
lems see lChow fc MonaghanI 1^93))- Despite the difficulties 
with modelling some boundary conditions, we believe that 
an entirely grid free method is the best solution to modelling 
these kinds of fiows. 

The advantage of SPH is that it handles shocks by pre- 
scribing an appropriate artificial dissipation rather than by 
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the use of exact Riemann solutions. Although such artificial 
dissipation can lead to greater dissipation than is the case 
for Riemann solver methods, the fact that the SPH resolu- 
tion adjusts for changes in density compensates because the 
dissipation terms axe proportional to the resolution. As a re- 
sult, one dimensional tests show that SPH could give results 
with similar acc uracy to either the Riema nn methods or the 
HLLE method dChow fc Monaghan|[l997il . In this paper we 
describe our algorithm in detail and concentrate our tests 
on shock tube problems in flat space. Application to both 
jets and disks and high energy ion collisions will be given 
elsewhere. 

Throughout the paper we use the Einstein summation 
convention, where Greek indices (^, v) indicate four vector 
components (0,1,2,3), and Latin signify spatial sum- 

mations only (1,2,3). We work in geometric units, where 
c = G = k = 1. And finally, we note that throughout our 
SPH equations we tag quantities associated with particles 
by the indices a and h. 



2 THE RELATIVISTIC GAS EQUATIONS 

We assume the gas consists of identical baryons with rest 
mass mo moving with four velocity U'^ in a space with met- 
ric coefficients g'^'^ . We can then completely specify the state 
of the fluid by defining this four velocity, the baryon num- 
ber density (n), the internal energy (e) and the isotropic 
pressure, (P). The energy momentum tensor can be written 



T^" = {n + ne + P)^^ + Pg" 



(1) 



These quantities are measured in the frame co-moving with 
the element of fluid. We assume that an equation of state 
specifying P in terms of the e and n is available together with 
a prescription for determining the composition for given en- 
ergy and density. In the present paper we assume the equa- 
tion of state is for an ideal gas of relativistic baryons how- 
ever, any equation of state can be applied. 

We are assuming here that the gravitational perturba- 
tion effects of the modelled ffuid can be ignored. This trans- 
lates to a background metric which is stationary. In order to 
decouple the time and spacelike integrations and interpola- 
tio ns, we use the 3-1-1 splitting of the space-time first done 
bv lArnowitt. Deser and Misnerl (Il962ll . This yields metrics 
of the form 

ds^ = ~{a - I3'fii)dt^ + iPidx'dt + rjijdx'dx^ (2) 

In this version of the equations, we further impose 
staticity upon the metric through enforcing that the shift 
vector (/3') is zero. The term a is the lapse function and is 
used to define the distance between the foliated, space-like 
hyp ersurf aces. 

We can then start by enforcing the conservation of par- 
ticles, 



(nt/");^ = 



(3) 



and that the divergence of the stress-energy tensor is also 
zero. 



rpfJ 



(4) 



From Equation Q 

dtN ^ -d,{Nv'), (5) 

where A'^ number density measured in the laboratory frame 
and is given by A'' = -yn and 

1 



7 



^ — 1 



(6) 



Equation 2] can be projected in directions parallel and 
orthogonal to the four- velocity (7^. The orthogonal compo- 
nents are given by, 

y'ihU^),, = .^{iNhU''U''g^,., - ftP) (7) 
iV7 z 



where 



h = l + e+ — 
n 



(8) 



and is the relativistic enthalpy. We then define the momen- 
tum density to be 



a 

=^«'7^(n + ne + P), 



(9) 



Here «* is the usual transport velocity dx' /dt, related to the 
four velocity through 



a 



(10) 



and y'^ is the determinate of the spatial part of the metric, 
and can be related to the full metric's determinant through 



The parallel component reduces to 



where 



dtE + dj{E-g°°P)v' 



E = 0^(n + ne + P)7^ - P. 



(11) 



(12) 



(13) 



The above equations are written in terms of momentum 
and energy per unit volume. In order to set up SPH equa- 
tions we need momentum and energy per baryon. However, 
due to the fact that we are in curved space, the spatial met- 
ric determinant y/rj needs to be incorporated into our spatial 
integrations. 

The momentum per particle q and the energy per par- 
ticle e then become 



and 



E 



VVN 



7(1. 



-e+^) 
n 



(14) 



(15) 



Using these definitions we find, using the Lagrangian 
time derivative 



— = dt + v^ dj , 
dt 

that the momentum and energy equations become 



N* 



(16) 



(17) 
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and 



dt iV 
with the continuity equation 

dN* _ 
dt ' 



(18) 



(19) 



In order to express these variable s, we have transformed to 
the number density N* {D* in ISieeler fc RiflFerd J2000l) '). 
which incorporates the volume element through A'^* = 

Making this transform alleviates the difficiulties expe- 
rienced by Laguna et al. (1993) with the approximations to 
spatial integrals in curved space. It allows us to use the 
flatspace kernels of traditional SPH (see Section 12.111 . It 
has also been shown that A''* satisfies the continuity equa- 
tions, whereas using results in extra terms being required 
JSieder fc Riffertil2000l) . Note also that there are no time 
derivatives of hydrodynamic variables on the right hand side, 
m aking the conservati o n equations more stable t han those 
of iLaguna etli] Jl993h iNorman fc Winkle jll983l . 

Particles are moved according to 



dt 



(20) 



These equations are identical in form to the non relativistic 
equations and they can be solved in the same way. 



2.1 SPH relativistic equations 

The baryon gas is replaced by a set of SPH particles with the 
number of baryons in SPH particle a denoted by Va- We use a 
and b as indices which should not be confused with vector or 
tensor quantities. Care is taken in the discretization process 
to ensure that the terms are symmetric in a and b, which 
not only increases calculation efficiency, but helps to ensure 
conservation of linear and angular momentum. 

Smoothed particle hydrodynamics is based upon the ap- 
proximation that 



< 



A{x) J yl(x')iy(lx'-x|,/i)dx'. 



(21) 



where A is any field variable, < A > is the approximation 
to A and W is some kernel, normalised such that 



J W{x.)dV = 1. 



(22) 



If we are in curved space, we can assume the same relations, 
except with the addition of y/rj, in that 



<Aix)>= / A(x')W'(|x' -x|,/i)0ydx' (23) 



and 



W(x)y^dV = 1. 



(24) 



If we then make the assumption that the W kernel can 
be described as some spatially varying function times the 
usual fiatspace kernel, W, 



W = /(x)VK 



(25) 



iLaguna et alJlTggSl) . then equation 12411 becomes 

j W(^)f(yL)^dV = 1. (26) 
which, if we compare to equation 1221 we can see that 



Equation (1231 then becomes 

W^(|x'-x!,/i) 



(27) 



< A{x) > 



A(x'; 



y/rjdyi 



N, 



(28) 



where i^i, is the number of baryons depicted by particle 6, 
and given by 

ub = ^NAV (29) 

Note how equation 1281 1 is identical in form to tradi- 
tional SPH, allowing us to use the usual SPH methods. 

The momentum, continuity and energy equations take 
the form: 

Conservation of Momentum: 



dt 




(30) 



Note: this equation is the same as that generated by 
iPrice fc MonaghanI ll200iri , although they used a variational 
principle. 

Conservation of Energy: 



dea 

-7r = -i^'"' 



P r,00,, P „00„ 

H r^To r ilab 



^,2 



and Conservation of Baryon Number: 
dN* V — • 



(31) 



(32) 



In the momentum and energy equations the quantities 
Ilab and Q.ab are dissipation terms which we discuss below. 
The equation of state used is that of the perfect fiuid 



(r - i)N€ 



(33) 



and is solved through a solving a no nlinear equation in 
the co nserved variables, as described in lChow fc MonaghanI 
1^23) ■ Here F is assumed to be constant, which allows the 
local sound speed, to be derived as 

h dn "n? de 
= ^-^^ (34) 

An issue which arises in the formulation of relativis- 
tic equations by using SPH particles is whether or not the 
kernels should take a Lorentz invariant form. However, it is 
clear that the representation of the gas by a set of SPH parti- 
cles is not an invariant process. For example, in our frame we 
might begin with the SPH particles arranged on a cubic grid 
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in analogy to the subdivision of space into cubical cells in a 
finite difference scheme. In another inertial frame the SPH 
particles will appear to be on a non-cubical grid, and the 
cells of the finite difference scheme will no longer be cubes. 
The computing observer in this other frame might decide to 
use either more or fewer particles and arrange them in some 
other way. The only issue is to calculate the dynamics as 
accurately as possible for given resources. 

Furthermore, we argue that the use of a flat space ker- 
nel is valid under the assumption that space is locally flat 
and that the radius of curvature of the background space 
time is large compared to the smoothing length. Should the 
scale of hydrodynamic variation (relative to the smoothing 
length) become small, particularly in an anisotropic fashion, 
then there is a clear a rgument for using spheroidal kernels 
iFulbright et al.lll995^ . The complications caused by this 
variation, in terms of relativistic consistency, are considered 
too great for this work and so are not attempted at this 
time. 



3 DISSIPATION TERMS 

A natural way to set up dissipation terms for SPH would 
be to use the continuum dissipation terms. In the case of 
non relativistic fluid dynamics it has turned out that more 
robust algorithms can be found which are not based directly 
on the contirmum dissipation terms (though they are equiv- 
alent when the continuum limit of the SPH equations is 
taken). In the case of relativistic fiuids the situation is more 
complicated because there are no satis factory dissipation 
terms t o begin with. T hos e proposed bvlLandau fc Lifshitd 
(|l95i), lEckard Jl940l) or IWeiriberd ^1972^ ar e know n to 
lead to i nstabilities (see 'I sraeli il97(j l. Ilsraell lll979ll and 
iHiscock fc Lindblom ( f985),). 

The dissipat ion terms were set up by 

IChow fc Monaehanl (Il997l^ in analogy to the dissipa- 
tion terms used in Riemann solver methods. These latter 
dissipation terms are based on adding to the average fiux 
from left and right cells a term of the form 



(35) 



where the A denote eigenvalues of the gas dynamic equa- 
tions, the e denote eigenvectors and r and i denote right and 
left states. These eigenvalues have the dimensions of veloc- 
ity and represent signal velocities. The F are the dependent 
variables of the differential equations: in the present case 
the relativistic momentum, energy and baryon number. In 
addition various limiters ma y be used to pre vent unwanted 
oscillations (see for example iLe Veau3 il997l) ). 

Following this idea the simplest dissipation term for the 
SPH momentum equation is 



Ilab = - 



KVsig{cia - qb) . j 



(36) 



where the eigenvalue is replaced by a signal velocity Vsig 
which we discuss below, and the eigenvectors are replaced 
by 

Tab 



Vab 



(37) 



which is the unit vector from particle 6 to a. 

The corresponding quantity for the energy equation be- 
comes 



Kvsig{ea — eb)j 



(38) 



In these expressions N*i, is an average of A'^^ and A^^ 
which in this paper we take to be (A^^ + N^)/2 and K is 
a constant. The quantity Vsig denotes the signal velocity 
which is the analogue of the eigenvalues used in the Riemann 
solver methods. We have also symmetrised the dissipation 
terms to guarantee that linear and angular momentum will 
be conserved. In this paper the dissipation terms are set to 
zero if the particles are moving away from one another. 

Although the form of the dissipation terms given above 
are consistent with the global conservation of linear and an- 
gular momentum and energy they are not consistent with 
the requirement that viscous dissipation should result in an 
increase in the thermal energy. Returning to and 1131 
and taking the special relativistic limit, we can write the 
thermal energy e in the form 



£ = 7(e- 

from which it follows that 
de de 

— = 7 ■ 

dt ' dt 



dq P d'y 
'dt N'dt' 



(39) 



(40) 



If s denotes the entropy per mass in the co-moving frame 
we can write 



yds de P dn 
dt dt n? dt ' 



(41) 



and substituting the non dissipative momentum and energy 
equations 1171 and 1181 together with 1391 the change in 
entropy is zero as expected. If we include the dissipation 
terms we find the rate of change of e for SPH particle a can 
be written 



deg 
dt 



TT-TtH 5— > rUbVab-VaWab 

N dt n2 ^ 



+ 7a 5^ fn-t, (HafcVa - Qab) ' V aWab 



(42) 



The first two terms on the right hand side are not associated 
with dissipation. The last term incorporates both viscous 
dissipation and thermal conduction. The rate of change of 
the entropy per mass is then 



='ya'^mb (HabVa - ilab) ■ ^ aWab 



dt 



(43) 



where the right hand side only involves the dissipation 
terms. In this expression for the entropy change s is mea- 
sured in the co-moving frame, while the time is measured in 
the laboratory frame. It is useful to write the rate of change 
of entropy per unit mass s as the four divergence of an en- 
tropy current 5''' = nsU^ according to 



d{Ns) 
'dt 



N— = - — = -V^ + V ■ (Nsv). 



(44) 



It is convenient to first consider the viscous dissipation 
terms. These can be isolated by setting the thermal terms 
to zero. We find from 1431 that the viscous dissipation term 
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{la - 76 - (7a Va " 76^6) ■ j (Va • j)) TabFab, 



(45) 



where we have replaced VaWab by VabFab where Fab ^ is 
a function symmetric in the coordinates of particles a and 
b. 

This expression is not positive definite because the 
terms from Uab involve velocities along the line joining the 
particles while 7 terms from the energy involve n^. To make 
the viscous dissipation positive definite we need an alterna- 
tive form of e in the dissipation terms. A simple and effective 
choice is the following. We first define v • j = «*, and 7* by 
l/^/l — (w*)^. We then replace in the dissipation term 
by e* where 



the viscous dissipation term then becomes 



Kvs 



mb- 



- V* 



i+Vb 



(46) 



rabFab, (47) 



To establish the sign of this expression we note with 



Vab 



we can form 

7ab = 



■Vg - Vb 

1 - vlvl 



= 7a7b(l - ^aVt) 



VI - (Kb)' 

and the dissipation term can then be written 

7a (1 - 7ab) rabFab- 



(48) 



(49) 



(50) 



Since 7*;, J5 1 and Fab ^ the viscous dissipation term is 
positive definite. The dissipation terms which we use in our 
equations are therefore Ilai,, which is unchanged, and 



Kvsig{e*a - e;)j 

Nab 



(51) 



The rate of change of entropy can be obtained from 
14411 when the only dissipation is due to viscosity and from 
the argument above, the change in the entropy from vis- 
cous dissipation is positive whether measured in the lab- 
oratory or the co-moving frame. Unfortunately there is 
no agreed form of the transformations of the thermody- 
namic quantities b etween dif f erent in ertial f r ames. Fur- 
thermore, alt hough 'Weinberd |l972), lEckartl lll940ll and 
iLandau fc Lifshitz ( 1959.) give entropy currents these de- 
pend on the forms they choose for the dissipation tensor. 
None o f the dissipation terms prop osed in the literature are 
stable jHiscock fc Lindblomlligia) . As a consequence it is 
not clear to us how to construct an appropriate entropy 
change from our dissipative terms when the thermal terms 
are included in the dissipation. 



3.1 Signal Speeds 

The si gnal speed Vsig was considered bv lChow fc Monaehanl 
^^23). In the non relativistic the speed of ap- 

proach of signals sent from SPH particles a and b towards 



each other iMonaghanll992t) . If these particles are separated 
by a distance r, and they are at rest, the signals meet at a 
time r/{ca + Cb) where Ca is the speed of sound of particle a. 
The signal speed is then taken as Ca + Cb- If the particles are 
moving terms must be added because the speed of sound is 
the propagation relative to the fluid which is itself moving. 
In the non relativistic case the signal speed becomes 



{Ca + Cb) - V, 



ab ■ J- 



(52) 



IChow fc MonaghanI (|^93) tested several generaliza- 
tions of this signal speed making allowance for the relativis- 
tic velocity addition formula. However, they were not able 
to derive a suitable Vaig. 

It is possible to be more precise about the form of Vsig by 
considering the sound waves sent between two SPH particles 
a and 6. The procedure we follow is to determine the wave 
number 4-vector of a sound wave from a particle a towards 
particle b in the co-moving frame of a. We then transform 
this back to our laboratory frame to get the apparent speed 
in the laboratory frame. We then repeat the procedure for 
particle b. Finally we construct Vsig by taking the average 
of these speeds. 

With this in mind, we take two inertial frames, the lab 
frame K, in which two particles a and 6 are moving, and 
the frame K' , in which particle a is at rest. At some time 
particle a emits a signal towards the second particle b. In 
K' , we know this signal to travel at the rest sound speed, 
Cs which can be calculated by thermodynamic quantities 
known in the co-moving frame. Our aim is to determine the 
apparent speed of sound Cs in the lab frame K. 

Let the wave number 4-vector of the emitted sound in 
the frame K he k'^ — {lu, k) where uj — kcs and k is di- 
rected from a towards b. In the frame K' of particle a the 
wave number 4-vector is (fc')'' = (aj',k') where ui' = k'cs- 
For convenience we rotate the laboratory axes so that the x 
axis of the lab frame is parallel to Va and denote the angle 
between the new x axis and the vector Yba from a to 6 by 6. 

The Lorentz transformation between the frames K and 
K' shows that 



7at^ 



;(0) 



(53) 



together with transformations of the wave number which 
lead to the aberration formula. From the invariance of fc^fc^ 
we deduce 



(54) 



where c is the speed of light. In the following we will assume 
that the speed of sound is scaled with the speed of light 
as in our other equations. Using the relation between the 
frequencies we can solve the previous equation for Cs. We 
find 



1(1 -c?)± 0.^(1 



1 — v^cii 



(55) 



where W|| is the component of the along the line 
joining the particles a and 6 and vx is the perpendicular 
component (all measured in the laboratory frame). One can 
immediately see that for the one dimensional case (where 
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Figure 1. Frames of Reference for soundspeeds: In the laboratory 
frame (on the left) we see two particles, a and b, travelling with 
unique velocities when a sends a signal to b. The right frame 
shows a at the origin, with the x-axis aligned with a's observed 
velocity in the lab frame. Particle a still transmits to particle b, 
but at a different speed. 



vx =0), this reduces to the usual relativistic addition for- 
mula 

Cs±V|| 



''I I 



1 ± CaV\l 



(56) 



Written in this way, equation I55II for the s ignal speed 
is the same as the eigenvalue A± deduced by iFont et all 
l)l994) for their acoustic waves travelling within a relativistic 
medium. By using this definition for our signal speeds (equa- 
tion H55|l 'l we ensure that the sound speeds remain causal 
and the artificial viscosity terms {Tlab and Qab) do not lead 
to acausal communication. 



4 TIMESTEP CONTROLS 

The algorithm used for these tests is based upon the 
predictor-corrector integrator. In order to show how this 
method works, we introduce 3 vectors. F is a five vector, 
containing the conserved quantities A''*, e and g', x is the 
usual location vector and H holds all the hydrodynamic vari- 
ables such as e, P, v and 7. We first predict F forward with 
an Euler step 

dF 

F = F + At^ (57) 

at old 

We then evolve the location vector forward through a 
second order prediction 



x = x + Atv+-(At)^- 



(58) 



For efficiency, we approximate the time derivative of the 
velocity with the known Tests were also run replacing 
equation 15811 with a standard predictor-corrector which pro- 
duced fundamentally the same results, indicating that this 
is a valid approximation. 

With these predicted values of F, we can solve for H 
using a Newton-Raphson root-finding method. The new H 
values allow us to generate new rates of change, These 
are used, in turn, to correct F through 



Finally, we can correct the values for H using the up- 
dated F. 

The timestep, dt, is given by 



dt = mm{Cvdtv, Cadta, Chdth} 
where C'a,Cv, Ch are constants and 



dta 



dth ~ min ■ 




(60) 

(61) 
(62) 
(63) 



The first two of these are the usual Courant conditions on 
the velocity an d the accelerat i on, re spectively. The third is 
that devised bv lThacker et al] i200Ci) and is a new condition, 
based upon the smoothing radius divided by the highest rel- 
ative velocity in that region. The algorithm seems rather 
insensitive to choice of the limiting constants and typical 
values chosen for Cv, Ca and Ch are 0.2, 0.35 and 0.35 re- 
spectively. 

It should also be noted that iThacker et al.l ll200(f) use 

the same 'Courant' conditions, but use 0.4, 0.25 and 0.2 as 
their limiting parameters C„, Ca, and Ch- The dominant 
limiting factor is usually dt^, and it has been further lim- 
ited in our scheme in an attempt to reduce overshooting 
due to the frequent supersonic motions and large gradients 
involved. 



5 NUMERICAL TESTS 
5.1 Boundary Conditions 

The hydrodynamic test problems presented in this paper re- 
quire a fixed box of known dimensions and so we employ a 
rectangular box. Two of the dimensions are periodic, and the 
third uses ghost-particles to maintain its structure. Period- 
icity is maintained by simply mapping the leftmost particles 
into pseudo-positions against the right boundary (and vice 
versa), and modifying the location vectors accordingly. 

The computational box used is defined as —0.5 ^ x < 
0.5 in the x-direction, and by the number of particles and 
the particle spacing in both the y and z directions. The 
initialisation procedure involves laying down the grid from 
X = —0.5 until close to the origin, where the particle spacing 
smoothly changes across to the lower density values in the 
right hand region. The upper x-boundary is then chosen to 
fit the particle spacing in exactly at around x — +0.5. 

The smoothing procedure involves smearing any discon- 
tinuity out across a few particle spacings in the x direction. 
To accomplish this, we simply define the particle position 
to be in one of three zones, the left steady state (Al), right 
steady state (Ar), or a transition region in between. The 
field variables A are then smoothed over the range by the 
rule 



A{x) 



Al + Arb^ 



(64) 



^ ^ At.dF dF , 
F = F H ( ) 

2 ^ dt dt old 



(59) 



where Aa; is the particle separation chosen to maintain 
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the relation between particles and desired density distribu- 
tion. The particles are all given the same baryon number i/i,- 
The initial (high density) particle spacing is chosen as an 
input parameter. This is then varied as the particles are ini- 
tially positioned, maintaining the relationship between Aa; 
and Vb to be the local baryo n number density (N*). T his 
is the same method used bv lChow fc Monaehanl il997ft in 
their 1-D application where it performed satisfactorily. 



6 SHOCK TUBE RESULTS 
6.1 Newtonian Shocks 

It is important that any new relativistic algorithm is first 
tested to ensure it gener ates the co rrect non-relativistic 
limit. The shock problem of lSodI il978h is probably the most 
recognised and so we have reproduced it here as a one dimen- 
sional problem, modelled in full 3-D. The initial configura- 
tion involves a discontinuity at a; = separating two states, 
the high energy left state {Ni = 1 x 10^, e; = 2.5 x lO"'^ and 
the right state {N* = 0.125 x 10 ^ = 2.0 x lO"'^) This con- 
figuration was initially done bv'Siegl er fc Rifferd i2000h and 
is replicated here, although calculated as a full 3-dimensional 
simulation. The numbers of particles in each coordinate 
direction (Numi) are given as ((Num^, = 2000, Num^ = 
8, Num^ = 8):(Num^ = 1000, Num^ = 4, Num^ = 4)) which 
corresponds to a computational box of 1.049 x 0.004 x 0.004. 
The smoothing lengths are chosen and evolved to maintain 
~ 57 neighbours. 

All of the following calculations are performed in the 
special relativistic limit, where a = 1 and 77'-' = rjij = 
diag{l,l,l}. 

Four graphs are shown, the number density, the 
isotropic pressure, the thermal energy and the velocity in 
the x-direction (the direction of shock propagation). Note 
the diffusion evident across the contact discontinuity in both 
Figure 13 and This can be reduced by increasing the res- 
olution of the simulation, which is currently limited by the 
constraints of single processor machines. In this particular 
example, the contact discontinuity is resolved over ~ 35 par- 
ticle spacings. These correspond to a computational time of 
~ 25. This time corresponds to ~ 26 sound crossing times 
which accounts for the scatter particularly evident in Figure 
|1| Figure|^shows some oscillation across the contact discon- 
tinuity, but no spiking as evidenced in other SPH discretiza- 
tions. The oscillation can be greatly reduced by tuning the 
viscosity parameter K, which in this case is 1.4, but this was 
not pursued to any great extent. 



6.2 Relativistic Shocks 

The next example we will show is the shock tube 
studie d in detail in one-d i mension bv [ Ha wley fit al] 
^1984^■ISchneider et alJ l|l993h . iMarti fc Mulleil l)l99(ih and 
ISiegler fc Riff erd (l200(t .A gain the higher density and pres- 
sure state is to the left, and lower energy states to the 
right. When the two states are brought together, a shock 
wave propagates into the low density gas to the right, and a 
rarefaction fan propagates leftwards into the higher density 
state. Between these two waves, there is a contact discon- 
tinuity where the pressure and velocity are continuous, but 




Figure 2. 3D simulation of the Newtonian Shock showing baryon 
number distribution (D) 




Figure 3. 3D simulation of the Newtonian Shock showing ther- 
mal energy, e 




Figure 4. 3D simulation of the Newtonian Shock showing veloc- 
ity in the x-direction 
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Figure 5. 3D simulation of tlie Newtonian Sliock showing 
isotropic pressure 



m - 




-0.4 -0.2 0.2 0.4 

Figure 6. ID simulation of the Relativistic Shock showing 
Baryon Number Density distribution 

the density drops. The initial conditions across the discon- 
tinuity are pressure; 13.3 : 1.0 x 10~® and Number Density; 
10:1. 

The baryon number density distribution is shown in 
Figure El where one can see that the locations of the shock 
front, contact discontinuity and the head and tail of the rar- 
efaction fan are well resolved. There is a slight over-estimate 
of the magnitude of the shocked shell. This is compared to 
Figure |7| which shows the identical 1-Dimensional phenom- 
ena, although the calculation is performed in 3 spatial di- 
mensions. Note that with the new Vsig and viscosity terms, 
there is no spike exhibited at the contact discontinuity. 

In moving from one to three dimensions there is signifi- 
cant increase in demand placed upon the dissipation terms. 
There is a lot more opportunity with the added degrees of 
freedom for the particles to interpenetrate, and for false com- 
munication to occur. Therefore we would not expect that 
what occurs in ID necessarily carries over into 3D. This is 
shown in the lack of resolution of the contact discontinutity 
and smearing of the shockfront. There is also a small amount 
of difTusion evident in the tail of the rarefaction fan. This 
is due to the artificial viscosity terms being unable to take 




-0.4 -0.2 0.2 0.4 

% displacement 

Figure 7. 3D simulation of the Relativistic Shock showing 
Baryon Number Density Distribution 
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Figure 8. 1-D baryon number density versus x location, with 
identical particle resolution in the x direction as Figure 171 

the directionality of the expansion into account. As a con- 
sequence, the expansion in x corresponds to particles also 
having y and z velocities. Because of the fixed boundaries 
these velocities correspond to colliding orbits in the y and 
z plane, which erroneously activate the artificial viscosity. 
The post-shock ringing is also seen in the dense shell. 

Many of these effects can be seen to be diminished by 
increasing the number of particles in the simulation. Due to 
the constraints of memory, we are currently unable to per- 
form the 3-dimensional calculation at the same x-resolution 
as the example shown in Figure |S] To highlight the issues of 
particle numbers, we then show Figure|H| which corresponds 
to a 1-dimensional calculation with identical x resolution as 
the Figure Q (350 particles in the x-direction across low 
density region) . Both Figure |7| and |S| show the simulation 
after 2000 timesteps. In the 1-dimensional calculation, the 
rarefaction fan is much better resolved, due to the viscos- 
ity terms not being activated. Notice also how there is no 
ringing in the shell region. The contact discontinuity is also 
better resolved, although less sharp than the high resolution 
Figure |S| which also displays a fiatter dense shell plateau. 

In order to examine the performance of the diffusion, we 
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-0,4 --0.2 0,2 0,4 

Figure 9. ID simulation of tlie Relativistic Sliock showing Pres- 
sure Distribution 




-U.I 1 : L u.i 

Figure 10. ID simulation of the Relativistic Shock showing Ther- 
mal Energy Distribution 

draw your attention to the 1-dimensional calculations shown 
in Figures 1^ and [T?71 and the three dimensional simulations 
shown in Figures [TTI and m\ 

In the one dimensional calculation, one can see the 
sharp resolution of the shock, and diffusion on the tail of 
the contact discontinuity. Also clear is the scatter in the 
pressure distribution at the same location. This spike is 
much curtailed in comparison to models which integrate the 
thermal energy equatio n rather than the total energy (see 
ISiegler fc Riffertl (l2000l) '). In 3 dimensions this spike IS re- 
duced to the same level as the surrounding scatter, although 
the diffusion caused by the large amounts of viscosity re- 
quired to stabilise the post-shock region is clear, highlighted 
by the excessive smearing of the contact discontinuity in 
Figure IT^ 

The curves most sensitive to scatter are the velocity 
curves, shown in Figure [T^ and 1141 (and IH . This is due to 
the manner in which the velocity components are calculated. 
This involves a rootfinder for the pressure term which then 
allows the magnitude of the velocity to be deduced. The 
magnitude of the velocity can then be decomposed into its 
constituent components through a comparison to the con- 
served quantity of momentum. In the one dimensional case. 



o : ^ ^ . I I 1 ^ ^ 
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Figure 11. 3D simulation of the Relativistic Shock showing Pres- 
sure Distribution 



-0,4 -0,2 0.2 0,4 

X displacement 



Figure 12. 3D simulation of the Relativistic Shock showing Ther- 
mal Energy Distribution 

this is a simple deconvolution, but is prone to more scatter 
in higher dimensions, as seen in Figure Fm 

The amount of scatter evident is due to the number of 
sound crossing times required to run the model. As pointed 
out previously, the computational box is defined by the num- 
ber of particles, not a physical distance. In this particular 
case, for this resolution this works out to be a box with di- 
mensions (1.05x0.02x0.02). If we look at Figure ITSl we can 
see that the sound speed is a high 0.72c in the high density 
region. For the time of 0.3865 units that the previous fig- 
ures are all produced at, we can see that this corresponds to 
nearly 15 sound crossing times. 

Continuing this analysis further is the Figure IT^ which 
shows the transverse velocities of the particles as a function 
of their x location. Note that these velocities are all under 
1.5 per cent of the local sound speed. 



7 SUMMARY 

In summary, we have shown here a set of relativistic SPH 
equations which correspond to a static background potential 
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Figure 13. ID simulation of the Relativistic Shock showing Ve- 
locity Distribution 



Figure 16. Velocities transverse to the direction of shock prop- 
agation 




X displacement 



Figure 14. 3D simulation of the Relativistic Shock showing Ve- 
locity Distribution 
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Figure 15. Sound Speed 



(space-time). They have no destabilizing time derivatives of 
hydrodynamic variables on the right hand side and have 
been shown to adequately capture basic shock flows. The 
shocks are captured through a new signal velocity term, 
which avoids the restricting requirement of an exact Rie- 
mann solution needed by some applications. The new deriva- 
tion of the artificial viscosity precludes any acausal commu- 
nication between particles and avoids the false spiking of 
other implementations. This false shock should be avoided 
as it may cause erroneous source terms in simulations where 
nuclear interactions, such as burning, occur. 

In their current guise the equations assume a stationary 
and static metric. The extension to metrics such as the Kerr 
is being studied. 
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